% Computes bi-variate transfer functions between two channels, given
% two input & one output channel
% Packs the residual abs(Eobs^2 - Epred^2) as the last two arrays
% calls r_resid to get the residual
% JUST COPY OF TRANSFN.M BU WITH SOME CHANGES TO
% ACCOMODATE ROBUST PROCESSING
% Latest date 12.03.2003

%zij = za zb - zc zd;
%tf = zxy zyx zxx zyy resEx resEy
function[tf] = r_trfn(SPM,ProcDef),

EQN=ProcDef.equations;
neq=ProcDef.neq;
A = size(SPM);
len = A(1);
tf = zeros([len neq]);

for l = 1:neq,
i = EQN(l,1);
j = EQN(l,2);
k = EQN(l,3);
p = sqrt(-1);

za = SPM(:,j,j);
if i>k,
   zb = SPM(:,i,k)+p*SPM(:,k,i);
elseif k>i,
   zb = SPM(:,k,i)-p*SPM(:,i,k);
end;

if k>j,
   zc = SPM(:,k,j)-p*SPM(:,j,k);
elseif k<j,
   zc = SPM(:,j,k)+p*SPM(:,k,j);
end;

if i>j,
   zd = SPM(:,i,j)+p*SPM(:,j,i);
elseif i<j,
   zd = SPM(:,j,i)-p*SPM(:,i,j);
end;

den = SPM(:,j,j).*SPM(:,k,k)- (SPM(:,k,j) - p*SPM(:,j,k)).*(SPM(:,k,j) + p*SPM(:,j,k));
top = [za.*zb-zc.*zd SPM(:,k,k).*zd-conj(zc).*zb];

tf(:,l) = (top(:,1)./den);
tf(:,l+2) = (top(:,2)./den);
end;

for i = 1:len,
data(:,:) = SPM(i,:,:);
tf(i,5) = r_resid(tf(i,3),tf(i,1),data,'Ex');
tf(i,6) = r_resid(tf(i,2),tf(i,4),data,'Ey');
end;
